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Till.  U.fclTli:  ThlOHY  OF  bXl’UXllON  INDUCE  WhimWG 
CdMl'UTLH  PROGRAM  GPIJCIFICATION 


fcy 

A.  N.  liieks,  .i.ISc.  and  J.  L.  McKeeman 

Thin  report  dene ribno  a  computer  program,  WORK  No.  6.1A,  which  evaluates  the 
* normu-l  mode*  equutionn  for  the  explosion  induced  whipping  motion  of  ships. 

The  first  part  of  the  program  determines  the  frequencies  and  shapes  of  the 
vertical  heaving,  pitching  and  whipping  modes  and  a  separate,  shortened 
version  of  the  program  (H.C.R.E.  No.  bill)  has  been  prepared  to  carry  out  just 
♦'ms  part  of  the  full  calculation.  The  programs  are  written  entirely  in 
FORTRAN  il  for  running  on  on  Atlas  computer.  A  version  is  also  to  be  written 
to  run  vita  the  1XJTRAN  i  compiler  on  a  KDF-9  computer.  As  the  proG^am  is 
rather  .long,  no  detailed  listing  is  given  here,  although  a  flow  diagram  has 
been  included  l Figure  l). 

J.  Tlio  derivation  of  the  equations  used  in  the  program  is  given  in  an 

earlier  report  [lj  and  the  equations  cover  only  the  simplest  possible  case; 
viz.,  I, he  completely  elastic  response  of  a  ship  not  too  close  to  a  pulsating 
but  non-migrating  explosion  bubble.  The  effects  of  migration,  i.e.  tne 
modified  magnitudes  of  the  later  pulses,  the  simple  change  in  the  explosion 
geometry  at  later  pulses  and  the  more  complicated  dipole  pressure  term 
induced  by  the  migration  ore  all  neglected.  Also  neglected  is  the  inter¬ 
action  of  the  spherically  diverging  flow  field  around  the  bubble  with  the 
finite  size  of  the  Bhip  structure.  This  effect  will  be  most  pronounced  for 
chorges  very  close  to  the  ship.  Multiple  images,  of  the  explosion  bubble 
due  to  the  proximity  of  the  water-surface  and  sea-bottom  can  however  be 
included  (up  to  a  maximum  of  four), 

3.  Although  these  complicating  effects  are  not  included  in  the  program, 
connideraole  provision  has  been  made  to  allow  for  their  later  inclusion  and 
the  program  os  it  stands  should  be  regJLrded  more  as  a  research  tool  to 
investigate  such  effects  than  as  an  established  means  for  predicting  whipping 
motions.  In  fact,  although  the  program  has  been  fairly  thoroughly  checked 
to  ensure  its  internal  consistency,  the  validity  of  the  equations  has  yet  to 
be  confirmed  by  comparison  of  the  predicted  results  with  experimental  data. 

The  purpose  of  this  report  is  to  present  a  clear  and  detailed  picture  of  the 
basic  program  before  the  proposed,  complicating,  modifications  are  carried 
out. 

k.  Unfortunately,  the  non-linear  nature  of  the  plastic  response  of  a  ship 
hinging  during  whipping  motions  precludes  the  possibility  of  using  a  normal 

/mode 
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mode  approach  in  such  cases.  The  more  direct  approach  outlined  in  [l]  would 

have  to  be  used. 

THK  EQUATION 

5.  To  make  the  report  as  self  contained  as  possible,  the  equations  are  given 
here  in  full,  although  with  little  explanation.  The  ship  is  considered  to  be 
representable  by  ’n’  lumped  masses  and  rotary  inertias  connected  by  (n  -  l) 
weightless  uniform  beams  of  equal  length  (whose  properties  vary  however  from 
beam  to  beam) .  The  elastic  and  buoyant  nature  of  the  ship  is  then 
represented  by  the  stiffness  matrix 


(A  +  K)  B 


C  +  K 


( 2n  x  2n) 


where  the  (n  x  n)  submatrices  A,  B  and  C  have  the  elements  defined  by 
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otherwise 


In  these  matrices,  l  is  the  length  of  each  beam  connecting  two  masses  and 
and  are  related  to  the  basic  beam  properties  by 


12(1  +  v)l. 
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L  and  v  are  Young's  Modulus  and  Poisson ’ a  Hatio  for  the  material  of  the  bean 

section.  I.  and  A.  are  the  moment  of  inertia  of  the  cross  section  and  the 
i  i 

cross-sectional  area  which  is  effective  in  shear  ( c.g .  for  a  solid 
rectangular  cross-section  a  x  o,  A  would  oe  P/3  ab).  The  stiffness  matrix  i> 
therefore  represents  fully  the  effect  of  shcur  deflections. 

6,  The  matrices  K  and  are  diagonal  matrices  representing  the  effect  of 
buoyancy.  It  is  assumed  tnat  when  the  ship  is  floating  in  its  equilibrium 
position,  the  sides,  at  the  waterline,  are  nearly  vertical.  'flie  diagonal 
element  k^  is  then  simply 

k  -  =  pA  .  ( i  =  1 ,  . . . ,  n ) 

i  vi  *  * 

where  A  .  is  the  waterplane  area  associated  with  the  ith  mass.  k  .  is 
wi  n 


given  by 


,  1  ,  .2  1  .  „2 

k  .  *  — —  k  .  £  »  —  oA  .  2, 
ri  12  i  12  wi 


( l  —  1^  * . ■ i  n) 


The  symmetric  matrix  S  is  transformed  into  the  symmetric  matrix  E  by  the 


transform 


E  =  M~2  S  M“2 


where  M~2  is  a  (2n  x  2n)  diagonal  matrix  with  elements  (}m)  given  by 


u-  *  (m.  +  m  . )”*  ) 

i  i  wi  ^ 


l  ) 


p.  =  (R.  +R  .)  2  ) 

n+i  a  wi  j 


(i  ®  I*  •  •  • ,  n) 


m.  is  the  mass  of  the  ith  lumped  mass*  m  .  its  added  water  mass  and  R.  its 

i  1  ’  Vi  i 

rotary  inertiao  R  ^  is  its  added  water  rotary  inertia  and  is  given  by 

1  2 
H  .  ■  ±7?  m  .  i 
wi  12  wi 

2 

7.  The  eigen-values  w^  (they  will  all  be  positive)  of  E  are  the  squares 

of  the  natural  angular  velocities  of  the  heaving,  pitching  and  whipping  modes 
of  the  ship.  The  frequency  of  the  ith  mode  is  then 

w . 

f.  *  ~ 

l  2v 

The  normal  mode  shapes  of  the  ship  are  related  to  the  eigen  vectors  of  E 
by  the  equation 


■  M’2  e. 

1  -i 


/Jiere 
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Here  y^  is  the  displacement  of  the  jth  mass  in  the  ith  mode  and  v^  is  the 

bending  rotation  at  the  jth  mass  due  to  the  ith  mode.  The  computer  program 
will  determine  and  print  out  the  nodal  frequencies  and  shapes,  the  shapes 
being  normalised  so  that 


*i  *i 

Mi  -1 

v.  v. 

-ij  J_-x_ 


(i.e.  I  {(m.  +  m.)y--2-*-(H-+R-)v.  .2}  =  1  ) 

j«l  1  wi ' “  ij  i  wi'  ij 

For  submarines,  k.  and  k  .  are  zero  for  all  i  and  the  first  two  natural 
*  i  rx 

frequencies  will  both  be  zero.  This  should  however  cause  no  difficulties 
since  the  eigen  value  routine  used  in  the  program  finds  accurately  orthogonal 
eigenvectors  even  for  repeated  eigen  values.  Physically  this  case 
corresponds  correctly  to  zero  pitching  and  heaving  frequencies. 

8.  The  normal  mode  equations  of  motion  are 

“i  +  wi2  °i  “  Ai  V  ^  ;  i  ■  1,  ...»  n  (1) 

In  this  equation  am  is  the  mode  coefficient,  VI t)  is  the  volume  of  the 
explosion  bubble  and  A^  the  mode  forcing  function  coefficient  defined  by 
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Mg  being  a  diagonal  matrix  with  elements  (“^)  given  by 


y.  ■  m  .  +  m  . 
x  wi  wx 


y.  -  R  .  +  R  . 
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) 
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where  m  .  ■  pA.l  ;  R  .  ■  m  .  I2  ;  A.  *  cross-sectional  area  of  the  water 

wx  x  *  wi  12  wx  ’x 

displaced  by  the  ship  at  the  ith  mass.  g.  and  h.  are  given  by 

J  J 
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where  the  geometries  are  indicated  in  Figures  2(a)  and  2(b).  For  a  given 
ship,  the  program  will  produce  a  series  of  sets  of  A..  corresponding  to  a 

given  series  of  attack  geometries. 

9.  For  non-migrating  bubbles,  the  volume-time  record  of  bubbles  can  be 
represented  fairly  accurately  by  a  single  universal,  non-dimensional, 
volume-time  record  v(x)  and  the  second  derivative  of  the  volume  can  be 
similarly  represented.  The  program,  at  present,  uses  linear  interpolation 

in  a  table  of  values  to  find  V(t).  Then  V(t)  ■  K  v(t) 


where  t  ■  t/T 

o 

For  a  charge  of  W  lbs.  of  T.N.T.  (or  the  T.N.T. equivalent  for  other 
explosives)  at  a  depth  D  ft.,  K  and  Tr  are  given  by 

K  «  1220  w1'3  (D  +  33)2/3  cu. ft. /sec. 2 

To  -  2.91*  V^/3/(D  ♦  33)5/6  secs. 

•  • 

The  present  tabic  for  v(t)  is  given  in  Table  1  and  is  shown  graphically  in 
Figure  3.  This  table,  punched  on  tape,  is  attached  to  the  end  of  the 
program  tape  and  may  easily  be  replaced  by  a  more  realistic  version  if  one  is 

available.  For  values  of  x  greater  than  the  maximum  given  in  the  table,  the 
program  puts  V(  t )  «  J. 


/10. 
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10.  As  commented  on  in  [lj ,  the  equations  of  motion  (l)  may  not  apply  in  the 
early  compressible  phase  of  an  explosion  and  V(t}  is  not  properly  defined  in 
this  phase  of  the  motion.  To  by-pass  such  difficulties,  the  compressible 
early  shock  motion  is  assumed  to  set  up  an  instantaneous  velocity 
distribution  along  the  ship,  and  the  function  V( t )  is  assumed  to  be  zero 
until  an  incompressible  flow  analysis  shows  that  it  would  otherwise  be 
changing  from  positive  to  negative.  The  first  1  alue  of  t  given  in  the  table 
corresponds  to  this  instant.  The  assumed  initial  velocity  distribution 
along  the  ship  is  given  by 


?c 

e 

V  . 

f " 

• 

V 

L'Cj 

-  M2 

h 

where  V  is  the  early  bubble  maximum  rate  of  volume  expansion,  assumed  to  be 
c 

achieved  instaneously.  From  an  incompressible  bubble  analysis  is  found 
to  be  given  by 


V 

c 


-  5950 


w^3 

(D  *  33)1/6 


cu.ft ./sec. 


The  program  will  accept  as  initial  conditions  for  the  integration  of 
equations  (l)  either  an  arbitrary  initial  velocity  distribution  upecified  on 
the  data  tape,  or  the  one  specified  above.  In  the  latter  case,  the 
distribution  will  also  be  printed  in  the  output.  With  either  form  of 
initial  velocity  distribution,  the  initial  conditions  for  the  equations  (1/ 
are 

a.  -  0  ) 

X  ) 

)  when  t  *  0 

) 

) 

For  the  standard  initial  velocity  distribution  these  are  equivalent  to 


6.  *  a. 
i  i 


*i 


v. 

'l 


V 

-c 


) 

)  t  -  0 

) 

) 

11.  Equations  (l)  are  integrated  in  the  program  by  a  step-by-step 
Runge-Kutta  integration  procedure  which  requires  the  selection  of  a  suitable 
tine  step  length  to  maint  in  accuracy  (or  even  stability).  For  stability 
it  is  necessary  that  the  incremental  time  step  should  not  exceed  about 
1/6  of  the  shortest  per: jd  (2w/w^)  inherent  in  equations  (l).  In  addition, 

the  forcing  function  Av(t)  is  sometimes  very  rapidly  changing  (e.g.  near 
each  bubble  pulse)  and  sometiaes  slowly  changing  (between  the  pulses). 


a.  *  0 
x 


0.  ■  a.  -  A.  V 
l  l  Ac 
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The  time  step  must  be  short  enough  to  represent  the  forcing  function 
accurately  even  when  it  is  changing  most  rapidly.  To  achieve  some 
economy  in  the  number  of  time  increments  used,  the  tine  increment  varies 
during  the  integration.  It  is  normally  set  to  6t  *  (l/Nth  of  the  shortest 
period  present  in  equations  (l)).  W  is  given  on  the  data  tape  and  for 
accuracy  it  should  normally  be  at  least  8.  If,  during  the  integration, 
the  program  finds  that  the  dimensionless  increment  between  the  current 
tabular  values  of  t  in  the  VTt)  table  is  less  than  6t/TQ,  then  it  reduces 

the  increment  to  this  smaller  value.  The  increment  is  increased  again  as 
soon  as  possible.  Thus,  between  the  pulses,  the  integration  step  is  tied 
to  the  stability  of  the  differential  equations,  and  near  the  pulses  it  is 
related  to  the  increments  necessary  to  define  the  forcing  function.  If 
the  standard  tabular  forcing  function  is  changed,  its  role  in  determining 
integration  step  lengths  near  pulses  should  be  borne  in  mind. 

12.  Since  many  more  incremental  integration  steps  will  be  made  than  are 
necessary  to  see  easily  the  form  of  the  results,  the  time  increments  at 
which  the  results  are  printed  is  independent  of  the  integration  step  length. 

A  print  increment  6t^  is  given  on  the  data  tape  and  the  program  prints  out 

t,  ou  and  am  at  the  nearest  integration  steps  to  the  print  intervals.  The 

output  is  therefore  not  quite  at  equal  time  intervals. 

13.  A  knowledge  of  the  mode  shapes  (jr^,  v^)  and  the  mode  coefficients  a. 

and  their  derivatives  ou  determine  the  whipping  motion  completely. 

Generally  however  it  is  more  convenient  to  know  the  total  motion  (displacement, 
velocity,  bending  moment,  shearing  force  and  deck  stress)  at.  specified 
points  along  the  ship.  This  motion  at  a  point  distance  x  from  the  bows  of 
the  ship  may  be  determined  from  the  equations 
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(2) 


/where 
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where  Ng  ia  the  number  ol‘  rnodnb  being  ouiih  iderod ,  1  it*  Hip  croan-ttnrt.  tonal 

moment  of  inertia  for  the  boom  section  at  the  point,  x  ami  y  is  t.lie  «l.\n tancr  of 
the  deck  from  the  neutral  axis  of  that  particular  beau  section.  Tin* 
quantities  y^(x),  M^(x)  and  LL{x)  arc  colled  the  'positional  mode  coefficient!*  * 

in  the  program  and  are  'no  dinplacenont ,  bending  moment  and  shearing  force  ut. 
the  point  x  viien  the  lumped  mass  displacements  and  rotations  are  defined  by 
the  ith  mode  shape  (y^,  v^). 

lh.  Since  each  beam  section  is  of  lenrrth  ,  x  can  ba  expressed  aa 


x  «  (j  -  i)t  x  where  0  $  "x  ft 

and  x  will  lie  in  the  jth  beam  with  shear  area  A j  and  moment  of  inertia 
1^.  Y^(x),  M^(x)  and  S^(x)  are  then  given  by  the  equations 
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where  y .  .  and  r- ■  have  the  same  meaning  an  earlier  (displacement  and  rotation 

A  J  Ij 

at  jth  mas 3  in  ith  mode)  and,  again  aa  earlier 

12(1  ♦  v)I. 


e  ■ 
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a.c 
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2EI . 

a  ■  —  ■  —  ■ 

!J(1  +  2c.) 

J 

The  sign  convention  for  these  quantities  is  shown  in  Figure  4. 

If  required,  the  program  will  form  and  print  out  both  the  positional  mode 
coefficients  and  for  each  print  time  step,  the  summations  defined  by  (2). 
These  quantities  can  be  compared  directly  with  experiment. 

ROTARY  IUKRTIAS 

15.  In  all  the  above  equations  it  has  been  assumed  that  rotary  inertia  plays 
an  important  part  in  determining  the  mode  shapes.  Where  it  is  thought  this 
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in  not.  the  t'tton  t  >ir>  order  of  fh«*  above  >4 1 1  m  t  toiin  |  genera  I  lv  .’it)  may  t»«* 
It*  1  vo, I . 

Thu  (.11  x  .’ii)  at.  iffiiena  matrix  N  111  replaced  Iiy  *  h«*  (11  *  11)  matrix 
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;;  -  ((a  ♦  k)  -  nr  n  ] 

Hie  (.’n  x  .’11)  matrix  i;  ia  replaced  lty  the  (11  x  11)  matrix 


•  •  e  i 

»■:  »  1  s  n  ’ 


*_  j  • 

where  Mj  1  in  the  (11  x  11)  diagonal  matrix  wit.lt  el  entente  given  by 

"  1 

Ii .  «  (in.  +  n^i  )  ( 1  -  1 . n) 


K  baa  only  11  eigen  va.lumi  amt  vectors. 

a 


From  each  eigenvector  e^,  the  displacements  of  the  lumped  man  a  on  in  the 
ith  mode  ahape  axe  determined  by 

* 

y  .  »  M,  e . 

-1  1  1 

This  equation  however  now  only  given  the  displacements.  The  rotation® 
must  be  determined  from  the  auxiliary  equation 

v .  -  -C  Vy. 

■I  ’*■  1 

The  normalisation  of  the  modes  in  thin  caae  is  such  that 


y-  m,  y.  -  1 
'1  1  *1 


n 

i.e.  .1 
J-l 


(m.  +  m  .  )y .  . 

1  wi  J 1J 


1 


Hie  mode  forcing  functions  axe  then  given  by 


A. 

1 


1  _y  *  1 

T^T  h  m2  6  "  TT7 


;■!  wi  wi  'ij 


where  in  the  obvious  (n  x  n)  contraction  of  M0.  Similarly  the  initial 

velocity  distribution  equation  is  a  simple  contraction  of  the  rotary  inertia 
case. 

lb.  The  program  incorporates  all  these  modifications  to  the  basic  equation 
when  rotary  inertias  are  to  be  suppressed. 
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IT*  The  program  contains  the  following  subvmit.  ins*  i- 

(a)  1'HIWT  1) 

H»ia  ju'lnts  out  mods  frequencies  and  shA|*n. 

(b)  HUNT  .'i 

Thia  routine  |>i  intn  ths  forcing  runotton  ooof f to vautn  whan  those 
ar«  i  squired. 

(o)  HK  (N)i 

Thia  carried  out  a  Rungs  Kuttn  integration  ott  n  not  of  * U *  first 
order  differential  equations,  it  requires  a  subroutine 
AUX  (1,  TA)  to  aat  up  valuaa  of  the  right  hand  aide  of  the 
equations. 

(d)  AUX  (I,  TA)  i 

This  routine  eeta  up  the  right  hand  aidoe  of  the  equations  solved 
by  the  routine  HK.  The  right  hand  sides  have  the  form 

•  a 

V  (t)  And  the  routine  n««d*  Another  routine  l»UCDK)l  (TA)  to 

givo  tha  current  value  of  V(t)  at  time  t  ■  TA.  When  the 
explosion  bubble  ie  migrating,  the  A^  values  will  also  be  time; 

dependant.  This  case  is  indicated  to  the  program  by  setting  D 
on  the  data  tape  to  minus  the  actual  charge  depth.  AUX  will 
then  call  in  an  additional,  routine  KORMD(TA)  which  will  determine 
the  current  depth  D(t)  of  the  bubble  and  then  set  corresponding 
A^  values.  In  tha  non-eugrating  case  the  same  A^  values  are 

retained  throughout  the  integration.  The  parameter  1  determines 
whether  the  current  entry  to  AUX  is  the  1st,  t'nu,  .ird  or  Uth  for 
the  integration  step.  On  the  2nd  and  hth  entries  AUX  adds 
|  DKLIAT  (DiKLTAT  is  the  current  integration  step  length)  to  TA. 

(e)  LAMDAj 

Qiwen  the  attack  geometry  d,  h,  W,  x^,  H  and  the  charge  depth  D, 

which  is  possibly  time  dependent  (set  by  FORMD  prior  to  entry  to 
LAMDA) ,  thia  routine  calculates  g  and  h,  from  these  and  the  mode 

shapes,  determines  thp  *  a lues  of  A^ ,  the  mode  forcing  function 

coefficients.  It  also  tests  n,  and,  if  this  is  negative,  i.t 
calls  in  PRINT  2  to  print  out  the  values  of  A^. 

(f)  MATDVIj 

This  routine  carries  out  the  matrix  division  required  in  the  evalua- 

s 

tion  of  S  ,  the  stiffness  matrix  in  the  non-rotary  inertia  case 

(g)  UCEIGN; 

This  routine  is  a  master  routine  for  calling  in  a  suite  of  sub¬ 
routines  which  carry  out  an  liGW  type  determination  of  the  eigen 
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valuoa  aiw!  vortm^  of  a  ayiMnfttrii'  matrix.  ’Hu*  routines  *r« 
vtoiu'iibod  in  ,)*•(.  ml  iu  [.’]  .  Only  tha  f  i  rat  aij;an  vaIura  (in 

ascending  order  of  frequency)  atr  found  «u»d  the  routines  will 
give  orthogonal  vectors  ovnii  for  repeated  .nol.a  (e„g.  tero 
frequency  pi  toll  orul  heave  for  a  submerged  submarine ) .  The 
suite  of  aubrouf. inan  coneint*  ofr- 

NCHT1)  -  Householder  Tridi  agonalie  At  ion 

I1CTD1U5  -  Tridi  a^oiiaI  Matrix  root  biaaction  routine 

STURM  -  Sturm  sequence  evaluation 

NCUET  -  Separation  routine  for  equal/nearly  equal  roota 
NCTH11  -  Invorae  Iterat  ion  procedure 

wCOHTll  -  Orthogonalieat.  ton  of  vectore  at  a  repeated  root 
NCltKTR  -  Hack  transformation  of  voctora 

(h)  FOHMD  (TA) ; 

This  routino  ia  called  in  by  AUX,  when  the  migrating  bubble  case  ia 
specified  on  the  data-tape  by  setting  D  to  -  (charge  depth). 

It  ehould  determine  the  depth  of  the  explosion  bubble  centre  for 
a  given  value  TA  of  the  time  aince  the  charge  detonated.  It 
should  then  call  in  the  subroutine  LAMDA  to  set  values  of  X^. 

At  present  the  subroutine  is  merely  a  dunxny,  leaving  the  value  of 
D  unchanged.  The  program  can  therefore  only  deal  with  the  non- 
migrating  cose. 

(i)  SECDER  (TA); 

a  e 

This  routine  determines  the  value  of  V(t),  the  bubble  volume 

acceleration  for  the  current  value  t  ■  i'A  of  the  time  since 

detonation.  It  also  determines  the  largest  time  interval 

«  * 

consistent  with  a  reasonably  accurate  definition  of  V(t)  where 

this  is  changing  rapidly  (near  the  bubble  pulses).  The  normal 

integration  interval  for  subroutine  RK  is  DEL.  DELW  is  the 

next  integration  step  which  will  be  used.  If  the  time  interval 
•  •  •  0 

defining  V  is  less  than  DEL  then  DELW  is  put  equal  to  the  V 
interval.  If  however  the  interval  is  greater  then  DEL  then 
DELW  is  set  equal  to  DEL.  The  present  SECDER  routine  confutes 

a  ’scale*  time  ^  to  reduce  TA  to  a  dimensionless  time  t.  It 

c 

interpolates  linearly  in  a  table  of  values  to  find  a  dimensionless 

acceleration  v(t)  and  converts  this  to  the  required  quantity 

•  • 

V(TA)  by  another  scale  factor  K.  T^  and  K  are  described  in 

•  • 

paragraph  8.  The  current  time  interval  needed  to  describe  V 

is  T  times  the  current  increment  in  the  table, 
o 
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16.  This  sin  vary  considerably ,  depending  on  the  output  required  of  the 
calculation.  Tabln  2  ahowa  a  schematic  baaio  data  tape,  together  with  the 
input  fonaata  and  unite.  Ilia  daaoription  in  Table  indicate#  where 
additions  or  subtractions  to  the  baaic  data  tape  ahuuld  be  made.  iiovoral  of 
tba  parameters  (e.g.  n,  M0 ,  etc.)  aerve  both  aa  baaic  data  and  aa  control 

paraasters.  If  ».et  to  aero  they  return  control  to  an  earlier  point  in  the 
program  jo  that  aeotiona  of  the  calculation  may  be  repeated  aeveral  time#  with 
different  sets  of  data.  The  program  can  only  be  properly  terminated  by 
using  these  control  parameters  to  return  control  to  the  start  of  the  data  tape 
and  than  setting  n  ■  0,  These  control  parameters  are  described  more  fully  in 
the  following  notes. 

IQCTS  OH  THE  DATA  TAPK 

19.  (a)  The  data  tape  ia  in  two  logically  distinct  parte:-  the  first,  from 

n  to  y  givaa  details  of  the  ship,  its  girder  strength  etc.,  and  the 
points  along  it  of  particular  interest.  The  second  port,  from 
°1  to  *F'  dives  details  of  the  attack,  geometry  and  the  integration 

procedure.  The  program  ia  arranged  so  that  the  second  part, 
from  ,  may  be  repeated  aa  often  a a  required  so  that,  at  one  run, 

several  different  attack  geometries  may  be  used  against  the  same 
ship  geometry ,  the  ship  normal  modes  being  evaluated  only  on 
the  first  calculation.  To  finish  calculations  for  the  one  ship, 
a  dummy  value  zero  should  replace  the  next  n^;  this  returns 

control  to  the  start  of  the  data  tape,  where  a  completely  new  ship 
can  be  considered  if  required  by  repeating  the  entire  data  tape 
with  the  new  ship  values.  Alternatively,  a  second  duntqy  value, 
zero,  for  the  new  n  will  stop  the  program. 

(b)  If  only  the  mode  shapes  and  frequencies  are  of  interest  and  these 
are  to  be  evaluated  for  several  different  ships  or  several  times 
for  the  same  ship  with  minor  changes  jin  the  ship  details,  this  is 
more  sisqply  accomplished  by  setting  Ng  *  0  in  each  set  of  ship 

data  and  omitting  the  rest  of  the  tape  for  each  case.  When  N2 

is  read,  if  it  is  zero,  control  is  returned  immediately  to  the 
start  of  the  tape.  A  separate  shortened  version  H.C.R.E,  Ho.  6lB 
of  the  program  has  been  made  which  evaluates  only  the  normal 
modes  and  frequencies.  The  data  tape  for  the  short  version  is  as 
for  the  full  program  down  to  but  then  repeats  if  more  than  one 
calculation  is  required. 
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(  c )  rpHitionul  tknto  Coefficients 

These  are  used  during  tao  integration  to  compute  the  total 
displacement ,  velocity,  bending  moment ,  shear  force  and  deck 
stress  due  tu  all  tao  included  modes  acting  together.  They  may 
however  be  of  interest  by  themselves  and  will  be  printed  out  if 
m  iu  specified  as  a  negative  value. 

(d)  'the  histories  of  some  of  the  quantities  displacement,  velocity, 

bending  moment  etc.  nay  not  be  required,  and  to  reduce  the 
output  only  the  relevant  ones  noed  bo  printed.  The  parameter# 
determining  the  print  are  ,  L0  and  I. 

■  1  displacement  and  velocity  are  printed 

■  0  "  M  "  w  not  printed 

L2  ■  1  bending  moment  and  shear  force  are  printed 

L2  -  0  "  11  ”  "  "  "  not  printed 

■  1  deck  stress  is  printed 

■  0  ”  M  ”  not  printed 

(e)  dumber  of  Charge  Images 

As  noted  in  paragraph  7t  the  free  water  surface  and  the  sea 

bottom  may  influence  the  ship  whipping  if  the  ship  or  the  charge 
is  close  to  either.  Their  effect  is  approximated  by  images; 
the  number  of  images  being  determined  by  | j • 

|n^j  ■  1  ;  charge  only 

In^l  *  2  ;  charge  +  surface  image 

|n^|  *  3  ;  ”  ♦  "  +  bottom  image 

In^j  «  4  ;  "  +  ”  +  "  "  +  image  of 

bottom  image  in  the  surface 

If  n^  is  set  to  zero ,  it  iranediately  returns  the  program  control 

to  the  start  of  the  data  tape  (n).  Although  these  images 
account  for  the  distortion  of  the  fluid  flow  due  to  the  sottom 
and  free  surface,  they  do  not  include  the  effect  on  the  bubble 
period.  Consequently,  for  charges  very  close  to  the  bottom  it 
would  be  better  to  put  |n^j  *  1  or  2,  ignoring  the  bottom  image, 

and  to  use  twice  the  given  charge  weight  for  W. 

( f )  Mode  Forcing  Function  Coefficients 

Apart  from  minor  changes  in  the  bubble  period  due  to  differing 
charge  depths,  the  extent  to  which  each  mode  is  excited  is,  by 
the  linearity  of  the  equations  (l),  proportional  to  the  mode 
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forcing  function  coefficients  Ju.  A  great  deal  may  therefore  be 

learnt  about  the  relative  severities  of  different  attack  Geometries 
simply  by  comparing  the  mode  coefficienta  for  the  geometries. 

Such  studies  minimise  the  number  of  complete  integrations  required. 
The  i.'e  will  be  printed  if  n^  is  given  a  negative  value.  To 

evaluate  several  sets  of  A^  for  different  attack  geometries,  the 

basic  geometry  set  d,  h,  W,  xw,  11,  D  may  be  repeated  as  often  as 
desired  if  the  values  of  V  are  all  made  negative.  jwj  will  then 
be  used  os  the  charge  weight.  If  any  positive  value  of  V  is 
encountered  the  program  will  continue  normally,  or,  if  a  zero 
value  is  specified  in  a  dummy  geometry  set,  control  will  be 
returned  immediately  to  the  data-tape  at  the  position  (see 
Note  (b)). 

(g)  If  D  in  the  attack  geometry  set  is  given  as  negative,  the  program 

will  assume  that  bubble  migration  is  taking  place  during  the 
integration  and  the  subroutine  FORMD  will  be  expected  to  give  the 
bubble  depth  at  any  specified  time.  This  will  allow  for  the 
change  in  the  position  of  the  bubble  between  pulses  but  will  not 
allow  for  the  effect  of  the  migration  on  the  magnitude  of  the 
pulses.  At  present  FORMD  is  a  dummy  routine  (see  subroutine  (h)). 

(h)  Initial  Velocity  Distribution 

This  should  normally  be  determined  by  the  equations  of  paragraph  9 
and  for  this  case  should  be  0  or  1.  If  ■  0,  the  initial 

values  of  a  and  §  only  will  be  computed,  as  initial  conditions 

for  the  Runge-Kutta  integration.  If  *  1,  the  actual  values  of 

£c  and  v^  will  also  be  computed  and  printed.  A  non-standard 

velocity  distribution  may  be  specified  on  the  data  tape  if 
required  by  setting  equal  to  2.  If  rotary  intertias  have 

been  included,  both  yc  and  vq  will  have  to  be  given  on  the  data 

tape,  but  if  rotary  inertia  has  been  neglected  throughout,  then 
only  yc  should  be  specified. 

OUTPUT 

20.  This  depends  very  much  on  the  data  tape  used  since  most  of  the  output  is 
optional.  The  output  has  been  extensively  titled  and,  apart  from  the  units 
used  is  self  explanatory.  The  maximum  output  consists  of 

(a)  Mode  shapes  and  frequencies.  The  displacements  are  in  feet  and 

the  frequencies  in  cycles/sec. 

(b)  Positional  mode  coefficients.  These  values  refer  to  the  ship 

under  unit  deformation  successively  in  each  mode.  The  units  are 
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displacement  -  ft.,  moment  -  lba.ft.,  shear  force  -  lbs.,  deck 
stress  (SIGMA)  -  lbs./sq.in.  This  last  quantity  is  the 
longitudinal  stress  at  the  specified  distance  'x*  along  the 
ship  and  height  y  above  the  neutral  axis. 

(c)  Attack  geometry  data.  Ihese  are  the  quantities  d,  h,  W,  j^,  H  and 

D  as  given  on  the  data  tape  except  that  negative  values  for 
W  and  D  (used  to  control  the  course  of  the  program)  will  have 
been  suppressed.  The  units  are  as  for  the  data  tape. 

(d)  Forcing  function  coefficients.  These  values  X^,  X2,  ...»  X^ 

are  in  (feet)  . 

a  • 

(e)  Initial  velocity  distribution.  y  is  in  ft. /sec.  and  v  in 


(f)  Time  histories.  For  each  print-time  increment,  the  time  (in  secs.) 
is  followed  by  the  values  of  the  mode  coefficients  a  and  §  (a). 
Then  for  each  specified  position  x,  the  displacement,  velocity 
etc.  due  to  all  the  specified  modes  together  is  given. 

MAXIMUM  SIZE  OF  CALCULATION 

21.  The  program  can  cope  with  a  ship  divided  into  up  to  kQ  lumped  masses 
and  kO  lumped  rotary  inertias  (i.e.  80  degrees  of  freedom). 

22.  Up  to  20  positions  may  be  specified  for  positional  mode  coefficients 
and  displacement  etc.  histories. 
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TABLE  1 

•  • 

DIMEUSIOMLESS  VOLUME  ACCELERATION  v( t ) 


0.16850 

+  0.0000 

2.18649 

-  2.0505 

0.19281 

-  0.5827 

2.28430 

-  1.6194 

0.28521 

-  1.9905 

2.34820 

-  1.1452 

0.31*729 

-  2.5879 

2.39739 

-  0.6150 

0.1*2805 

-  3.1353 

2.43765 

-  0.0106 

0.55009 

-  3.6433 

2.47166 

+  0.6945 

0.716U 

-  3.9090 

2.50104 

+  1.5392 

0.76462 

-  3.9090 

2.52688 

+  2.5837 

0.93061* 

-  3.6433 

2.55021 

+  3.9240 

1.05268 

-  3.1353 

2.57237 

+  5.7217 

1.1331*1* 

-  2.5879 

2.56638 

+  7.1310 

1.19552 

-  1.9905 

2.59759 

+  8.2644 

1.28792 

-  0.5827 

2.60543 

+  8.9034 

1.31223 

+  0.0000 

2.61352 

+  9.2785 

1.37936 

+  2.5207 

1.1*1973 

+  6.0896 

2.62008 

+10.4291 

1.1*1*802 

+12.9140 

2.62490 

+  9.8540 

1.1*581*6 

+19.2760 

2.63518 

+  7.7911 

1.46971 

+37.2590 

2.65371 

+  4.6281 

1.47366 

+52.6290 

2.66761 

+  3.1751 

1.47769 

+79.1790 

2.68312 

+  2.0977 

1.47859 

+85.4500 

2.72113 

+  0.5882 

1.47930 

+89.5690 

2.77276 

-  0.4515 

1.47962 

+91.0060 

2.80650 

-  0.8710 

1.48006 

+92.4780 

2.84889 

-  1.2455 

1.48036 

♦93.0770 

2.90710 

-  1.5855 

1.46449 

+  9.2785 

1.49258 

♦  8.9034 

1.50042 

♦  8.2644 

1.51163 

♦  7.1310 

1.52564 

+  5.7217 

1.54780 

♦  3.9240 

1.57H3 

♦  2.5837 

1.59697 

+  1.5392 

1.62635 

♦  O.6945 

1.66036 

-  0.0106 

1.70062 

-  0.6150 

1.74981 

-  1.1452 

1.81367 

-  1.6194 

1.91152 

-  2.0505 

2.03656 

-  2.2490 

2.06145 

-  2.2490  j 

tion  etepe  per  cycle  of 
ai(be«i  frequency  Bode. 
(*««  nere^repb  10) 
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FLOW  DIAGRAM 
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COEFFICIENTS 


INTEGRATION  OF  THE  DIFFERENTIAL  EQUATIONS 


